The H-mode has been considered as part of the standard operating scenario in the International Thermonuclear Experimental Reactor (ITER) due to its better energy confinement.^[1] Strong gradient in H-mode plasma routinely gives rise to the edge localized modes (ELMs), which eject a great number of particles with high energy onto divertor targets, and thus resulting in possible erosion of the divertor plates.^[2] Great effort has been made to understand the physics and better control of the ELMs via various external interventions, such as pellet injection (PI),^[3] resonant magnetic perturbation (RMP),^[4] impurity seeding,^[5] etc. Specifically, supersonic molecular beam injection (SMBI) has been a promising technique for ELM mitigation.^[6] The deposition location^[7] and the cooling effect^[8] of SMBI are relevant to ELM mitigation. According to the experimental results on HL-2A, it is found that the performance of ELM mitigation by the helium (He) is better than that by the deuterium due to the fact that the duration of ELM mitigation with helium SMBI (∼ 40 ms) is much longer than with deuterium SMBI (∼ 10 ms).^[8]

Gas injection is a good technique for externally interfering in tokamak plasma behaviors. High-pressure helium gas injection has been used for mitigating the disruption on EAST tokamak.^[9] Particle fueling can provide the dominant control of the size of the H-mode density barrier.^[10] Helium gas puff was utilized to achieve a two-dimensional (2D) measurement of edge turbulence with high temporal and spatial resolution, this diagnostic technology is named the gas puff imaging (GPI).^[11] The comparison of experimental GPI emission data with the simulation results can be used to validate the neutral transport code. The behavior of neutral helium in gas puff imaging experiments has been investigated.^[12] A coupled plasma-neutral fluid model has been used to study the interaction between the divertor plasma and the hydrogen gas target, and the simulation results agree well with the experimental data,^[13] suggesting that the fluid method can deal with the neutral transport in tokamak.

Kinetic simulation is needed for neutral transport when the fluid method is invalid. The solution of the Boltzmann kinetic equation for neutrals can be obtained via different methods. The direct simulation Monte Carlo method (DSMC) approach^[14] statistically mimics the behavior of real neutrals by tracking a group of model particles. Neutral transport can also be given by solving the Boltzmann equation with a non-linear BGK approximation.^[15] Currently, several well-known code packages are often used to study the plasma and neutral gas dynamics in the non-core region of a tokamak, such as EMC3-EIRENE,^[16] EDGE2D-EIRENE,^[17] B2-EIRENE (SOLPS),^[18] SOLPS-ITER,^[19] UEDGE-DEGAS 2.^[20] The continuum modeling techniques such as the Navier–Stokes equations are inadequate if the neutral gas is sufficiently rarefied. However, the rough fluid method is still applicable to neutral gas jet simulation^[21,22] in the edge high-collision region. The BOUT/BOUT ++^[23–25] code is a framework for the 2D and three-dimensional (3D) fluid simulation in curvilinear geometry. Many physical modules have been developed, such as the trans-neut module^[26,27] for studying macroscopic transport dynamics in 3D simulations. Recently, a test particle module is developed under the BOUT ++ frame work for studying the impurity migration patterns.^[28] There are relatively few studies devoted to quantitative information about the transport characteristics of helium species injected by the SMBI and the relevant plasma response. For simplicity, a fluid model is adopted for both plasma and neutral in this paper, and implemented in the BOUT ++ framework.

The transport of helium species in the edge plasma is investigated by using the newly-developed trans-imp code in the frame of 3D plasma fluid simulation code BOUT ++ in the present paper. The effect of the helium species on the edge-tokamak plasma is also studied by using the trans-imp module. The complete transport dynamics for the plasma and impurity components are divided into three different stages: steady phase before the SMBI, injecting phase during the SMBI and relaxing phase after the SMBI. The physical model, including transport equations, the initial and boundary conditions, is described in Section 2. Simulation results are illustrated in Section 3. Finally, the primary results are summarized in Section 4.

A simplified model is developed for studying the fundamental transport properties of the gaseous impurity SMBI into the deuterium tokamak plasma by assuming the quasi-neutral condition , where n is the charge number and Z is the atomic number. For simplicity, both the evolutions of the electric and magnetic field are ignored. All the drift velocities are omitted, and only the parallel fluid velocities are considered in the present model.

2.1. Transport equations

The plasma evolutions are determined by the following equations:

where N_i, V_∥i, and T_i are the number density, parallel velocity, and temperature of the deuterium ion, T_e is the electron temperature, is the parallel classical ion viscosity, the parallel classical ion thermal conductivity, the parallel classical electron thermal conductivity, , , and are the perpendicular diffusion coefficients of ion density, deuterium ion temperature, and electron temperature, respectively, subscripts Z₀ and Z_n denote the neutral particles and impurity ions, respectively, and P_j = k_BN_jT_j is the pressure for the j species (j = e, i, Z₀, Z_n) with k_B being the Boltzmann’s constant. The frictional force is written as follows:^[29]

The evolutions for the impurity are determined by

where N_Z0 is the impurity atom number density, V_Z0x is the impurity atom radial velocity, N_Zn, V_∥Zn, and T_Zn are the number density, parallel velocity, and temperature for the n-th (n = 1,2,…n_max, n_max = 2 for helium) charged impurity ions, respectively. The symbol ∇_x refers to the radial component of the gradient, and M to the particle mass. For all the He ions, the diffusivity values are assumed to be the same. Both the particle diffusion coefficient and the thermal diffusion coefficient are set to be 1.00 m²/s in the present work. The transport of the impurity ions in multiple ionization states within the non-uniform plasma is governed by several simultaneous processes, such as diffusion, convection, ionization, and recombination of the impurity ions.

In the present study only the dominant inelastic collisional reactions, including ionization and recombination, are taken into consideration. The momentum and energy source terms caused by inelastic collisions for particles with ionized state n are written as follows:

where subscript δ denotes inelastic collision, such as ionization, recombination, and charge exchange, refers to the particle number source term caused by the δ-type reaction that the particles in the n-th (here, n can be zero) order ionization state take part in, p_n = M_nV_n to the fluid momentum to the fluid energy. These cross-sections or rate coefficients are extracted from the atomic data and analysis structure (ADAS). For simplicity, the particle species symbol is omitted. The source terms are given by the following expressions:

where ∑_δ is the summation overall inelastic reactions. There are four possible values for δ(+1, −1, d, and r). The first two values ± 1 indicate that particles participate in the ionization and recombination reactions, respectively. The last two values d and r indicate that particles participate in charge exchange reactions as donors and receivers, respectively. I(n,δ) = 1 or −1 is a weighting factor to show a positive or negative contribution from a reaction, respectively. The relationship between the heat source S^h and the energy source S^ε is

We assume that the plasma does not absorb its own radiation. The total radiation power loss originates from the ionization, recombination, bremsstrahlung, and cyclotron radiation and is given below

Four kinds of corresponding radiation energy loss densities (all in units of W·m⁻³) are included in this paper and are shown below:^[30,31]

where χ_n is the ionization energy in unit eV, B is the magnetic induction in unit Tesla, and α_n and β_n are the rate coefficients (in unit m³·s⁻¹) for ionization and recombination, respectively.

2.2. Initialization and boundary conditions

In order to validate the reliability of the present model, the magnetic configuration (as shown in Fig. 1) with a circular cross-section is used. The initial H-mode profiles of the deuterium ion density, electron temperature, and deuterium ion temperature (shown in Fig. 2) are given as follows:

where x = (ψ_ped − ψ_N)/Δψ, with ψ_N being the normalized magnetic flux, ψ_ped, and Δψ being adjustable input parameters; a and b are also the adjustable input parameters; tanh is a hyperbolic tangent function. The location and width of the pedestal are set to be ψ_ped = 0.9 and Δψ = 0.025. In the present simulation, the parameters are set to be a = 3.8, b = 0.015, N_xin = 0.30 N₀ (where N₀ = 1 × 10¹⁹ m⁻³), N_{x out} = 0.01, N₀, T_{i(e),x in} = 250 eV, and T_{i(e),x out} = 1.0 eV. The radially dependent diffusion coefficients , , and , are calculated to keep the H-mode profiles stable before SMBI. When calculating the particle diffusion coefficient, imposed are the conditions of constant fluxes that the inflowing fluxes at core Γ_in and the outgoing fluxes at edge Γ_out are the same as those anywhere in the perpendicular direction. Under those conditions of constant fluxes, it is easy to keep the initial profiles stable. For example, Γ_in=Γ_out = Γ_⊥ = −D_⊥∇_⊥ N_i are set for keeping profiles stable, where Γ_in = −D_{⊥ in}(∇_⊥N_i)_in is a constant at core. Thus, D_⊥ = −Γ_in/∇_⊥N_i, while ∇_⊥N_i is calculated according to the initial density profile. The perpendicular thermal diffusion coefficients and are calculated by the same method. Strictly speaking, particle and energy transport in tokamaks cannot be described in terms of diffusion alone. A critical physical model with turbulent diffusion and the particle pinch, are used to explain the spatiotemporal dynamics of anomalous particle convection reversal.^[32] But the particle pinch is not considered in this paper due to the difficulty in dealing with the particle pinch. The edge transport barrier is characterized by these diffusion coefficients at ψ_N = 0.9 as shown in Fig. 2.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Contour plot of equilibrium poloidal magnetic flux.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (a) Initial profiles of plasma temperature and thermal diffusion coefficient, and (b) initial profiles of plasma density and particle diffusion coefficient.

On the inner boundary, the flux-driven boundary condition is imposed for plasma density(temperature). Neumann boundary conditions are set for all other evolving quantities. At the outer boundary, Dirichlet boundary conditions are set to be N_i = 0.01 N₀, T_i = T_e = 1 eV for plasma density and temperature, and Neumann boundary conditions are set for all other evolving quantities except N_Z0 and V_Z0x. For modeling the SMBI, the impurity density, and radial velocity are expressed respectively as

where Δθ = l_θ/r₀ = θ₁ − θ₀, θ_m = (θ₁ + θ₀)/2, r₀ and l_θ are the minor radius of plasma and the poloidal length of the supersonic molecular beam injected into the outermost magnetic surface. The parameters of He SMBI are set to be , m/s, Δθ = 0.628 rad, θ_m = 3.14 rad, and Δt = 0.3 ms (pulse width). The injection channel is located at θ_m = 3.14 rad. The shortest pulse width of SMBI is 0.2 ms–0.3 ms.

The interaction between He SMB and MHD turbulence in tokamak plasma plays an important role in fueling tokamak. The present study focuses on the fundamental transport properties for the background plasma and helium species. Figure 3 shows that the plasma is in the steady state before the SMBI. In the simulation, the He atoms are injected into the steady-state H-mode plasma through the SMBI. The spatiotemporal evolution of the plasma and impurity profile during SMBI are simulated by the newly developed trans-imp code.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Time evolution of (a) background plasma density and (b) temperature before He SMBI at ψN = 0.9 and θ = 3.14 rad.

All electrons are stripped off if the helium species are injected into the high-temperature region (see the yellow curve in Fig. 4). The accumulated He²⁺ ions become the main helium species with the He atoms and He⁺ ions decreasing. The SMBI is turned off after ∼ 0.3 ms. Without the particles leaving from the innermost or the outermost boundary flux surface, the total number of impurity particles keeps constant. The He²⁺ ions number reaches a peak value of 1.1 × 10¹⁹. Because the recombination reaction rate is much lower than other atomic and molecular reaction rates in the hot plasma (T_e ≫ 1 eV), it is difficult for He²⁺ ions to recombine into He⁺ ions. In the whole simulation, the total particle number of deuterium ions amounts to 6.0 × 10¹⁹, and keeps constant (see the purple curve in Fig. 4).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (a) Time evolution of total particle number of He atoms (), He+ ions (), He2+ ions (), He species (); (b) time evolution of total particle number of deuterium ions , with ∫Ω dV denoting volume integral and Ω being region to be solved.

3.1. Helium transport during SMBI

Particle, momentum, and heat transport of He species injected into a steady-state H-mode (t = 0 ms) plasma in both radial and poloidal directions are studied. When the SMBI is turned on, the He SMB front begins to propagate inward and towards the pedestal, then propagates outward, and finally stays at a position near the middle of the pedestal (ψ_N = 0.9), termed as deposition location, which is shown in Fig. 5(a). Along the He SMBI path, the local plasma is cooled down due to He ionization and radiation. The injecting He atoms are highly ionized at the propagating front where they have the highest ionization rate, while they have less ionization rate along the injection path. It leads the front layer of He SMBI to be highly ionized as observed in the experiment. The deposition of He SMBI in this paper is located at ψ_N = 0.92. The electron temperature at the deposition location (about 100 eV) is larger than the ionization energy of helium particles (e.g., the first-order ionization energy 24.6 eV and the second-order ionization energy 54.4 eV). It is inferred that He atoms are found in the whole injection channel during the early stage of He SMBI (see Fig. 5(b)). Subsequently, the radial profiles of He atoms and He⁺ ions tend to be more peaked, which means that the He atoms and He⁺ ions are completely ionized into He²⁺ outside the deposition location due to the high electron temperature and density. The He⁺ number density N_Z1 at the deposition location decreases with time going by, while the He²⁺ number density N_Z2 increases rapidly with time elapsing. Figures 5(d) and 5(e) show that it takes about 0.05 ms for He²⁺ ions to become as hot as deuterium ions. However, the temperature of He⁺ ions will become similar to the second ionization energy of helium after t = 0.05 ms, because He⁺ ions are lost through its ionization before reaching a high temperature of more than 54.4 eV. The pressure of He⁺ ions will not turn too high (see Fig. 6(a)), while the pressure of He²⁺ ions becomes more and more peaked inside ψ_N = 0.9 with time increasing (see Fig. 6(b)). The ionization of helium increases the density of electrons, He⁺, and He²⁺ ions. Initially, the density of He⁺ ions is very low, therefore, the first-order ionization is more intense than the second-order ionization. Then the increases of electron density and He⁺ ions density, but the little change of helium atom density, make the second-order ionization stronger than the first-order ionization. As a result, the pressure increases at an earlier time and then decreases. Note that the density peak of He²⁺ ions is located outside ψ_N = 0.9. The physical picture of the interaction between He SMBI and background plasma is that the He particles run inward because of inertia and collide with the hot electrons to produce He⁺ ions, continuously complementing the He⁺ ionization loss.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Time evolution of radial profile of NZ0 (a), NZ1 (b), NZ2 (c), TZ1 (d), and TZ2 (e) along injection path at the fixed poloidal angle during He SMBI.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Time evolution of radial profile of PZ1 (a) and PZ2 (b) along injection path at fixed poloidal angle during He SMBI. The normalized pressure unit is set to be P0 = 16.02 Pa.

Then the He⁺ ions further collide with the hot electrons to produce He²⁺ ions. The newly produced He cations exchange momentum with the background plasma and absorb energy from plasma.

The background plasma response to the helium SMBI into the H-mode pedestal is also investigated. The result indicates that the electron temperature and deuterium ion temperature decrease slightly with the increase of cold He particles injected by SMBI, but the deuterium ion density decreases sharply at the deposition location as shown in Fig. 7. Before the He SMBI, both the electron and ion poloidal temperature profiles are the same and uniform. After the He SMBI, electron temperature drops faster than the deuterium ion temperature at the deposition location (see Fig. 8). This is mainly because of the dilution effect of newly generated cold electrons and the absorption of the He impurity ionization reactions. However, the poloidal temperature profile of electron and deuterium ion exhibit different levels of flatness, which indicates that for electrons, their parallel thermal conductivity contributes more to heat transport than the thermal convection, contrary to which is the result for deuterium ions. The poloidal profile of deuterium ion density drops locally in the region of He SMBI, but peaks at the opposite position of He SMBI. This is because of the opposite velocity directions resulting from the frictional force between the main ion and the impurity ion. The impurity ion velocities are driven by the poloidal pressure gradient.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Time evolution of radial profile of Te (a), Ti (b), Ni (c) along injection path at fixed poloidal angle during He SMBI.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Time evolution of poloidal profile of Te (a), Ti (b), log10Ni (c), and V∥ i (d) at deposition location during He SMBI. θ is normalized by 2π.

3.2. Helium transport after SMBI

The remaining atoms and He⁺ ions will be quickly ionized in about 0.05 ms after the He SMBI. The radially localized peaked impurity ion density (N_Z2 ∼ 1.7N₀) can be several times larger than the local deuterium ion density (N_i ∼ 0.4N₀), and spread out of the deposition location in both radial and poloidal directions. The result suggests that the He⁺ ion temperature T_Z1 almost equals the deuterium ion temperature (see Fig. 9(d)). It is difficult to understand that the helium particles can withstand so high temperature without being completely ionized, which may be caused by the initial assumption that the impurity density cannot be lower than 10⁻¹⁰N₀ in the present simulation. In fact, N_Z1 near the core plasma should be close to zero. It takes about 0.05 ms for He²⁺ ion temperature T_Z2 to reach a steady value in the radial direction after the SMBI and to become as hot as deuterium ions because of the intensive energy exchange (see Fig. 9(e)).

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Time evolution of radial profile of NZ0 (a), NZ1 (b), NZ2 (c), TZ1 (d), and TZ2 (e) along injection path at fixed poloidal angle after the He SMBI.

The mean main plasma profiles return to the initial states in the radial direction after the He SMBI at t = 0.8 ms (see Fig. 10). During the relaxing phase at the deposition location, the electron temperature is lower than the deuterium ion temperature due to the newly-born cold electrons and the radiation loss of the electrons. After some time, the electron temperature becomes close to the deuterium ion temperature (see Figs. 11(a) and 11(b)). It takes about one millisecond for the deuterium ion density to become poloidally uniform via poloidal transport and the temperature will become uniform at t = 1.3 ms (see Fig. 11).

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. Time evolution of radial profile of Te (a), Ti (b), and Ni (c) along injection path after He SMBI at the fixed poloidal angle.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. Time evolution of poloidal profile of Te (a), Ti (b), log10Ni (c), and V∥ i (d) at deposition location after He SMBI. θ is normalized by 2π.

A simplified model is set up by reducing the moment equations of multi-component plasma. This model is utilized to investigate the transport properties of helium species and the response of the background plasma in the tokamak H-mode edge plasma when the He SMBI is used. In the model, the ionization and recombination of He species are taken into account but the charge exchange reaction is not considered. The model is implemented by a new physical module, trans-imp, in the frame of BOUT ++. The magnetic field configuration with a circular cross-section is chosen to simulate the helium species transport during the He SMBI. The results indicate that the impurity particles peak at the plasma edge at the beginning of He SMBI, existing in the form of completely stripped helium ions. In the injecting phase, the He²⁺ ions will be locally deposited at ψ_N = 0.9 where the electron temperature is higher than 54.4 eV when the source terms become more dominant than the transport terms. During the relaxing phase, the impurity ions are pushed into both inner and outer region by the large density gradient. In the future, our work will be extended in the following aspects: (i) further in-depth study of He-SMBI-plasma in diverter configuration, (ii) noble gas injection besides helium injection, and (iii) ELM mitigation after the He SMBI.